clear all
syms x y
w = 1+3i;
w = w/abs(w)*.9;
lambda =-1+2i;
lambda0 =  lambda;

q12_re = 0:1:10;
q12_im = q12_re;
Q = [];

for a = 1:length(q12_re)
    for b = 1:length(q12_re)
        Q(end+1) = q12_re(a) + j*q12_im(b);
    end
end
X = {}; Y = {};    
for i=1:length(Q);
    q12 = Q(i);
    q = (1-abs(w)^2)/(2*abs(lambda)^2*abs(w)^2)*abs(q12)^2 - real(lambda^2*w*q12)/abs(lambda)^2/abs(w)^2;
    r = solve(-x^2+y^2 +real(lambda^2)+q+real(w*q12), ...
        -(x^2+y^2)^2 +abs(lambda)^4 +2*abs(lambda)^2*q +2*real(lambda^2*w*q12) +(1-abs(w)^2)*(q^2-abs(q12)^2),x,y);
    X{end+1} = r.x;
    Y{end+1} = r.y;
end
XX = [];
YY = [];
for i = 1:length(Q)
    XX(1:8,i) = double(X{i});
    YY(1:8,i) = double(Y{i});
end
%%
for f =1:8
    figure(f);
    hold on
    for i=1:length(Q)

        if (imag(XX(f,i))==0) && ( imag(YY(f,i))==0)
            plot(XX(f,i),YY(f,i),'+');
        end
    end
end